function r2 = model_r2_mean_1D_matrix_form(para,t,parameters)
% para(1): nu
% para(2): delta^2/nu
% output: roughness square

mode = parameters(1);
L    = parameters(2);
a    = parameters(3);
freq = parameters(4);
dx   = parameters(5); % Aggregation length, unit site
l = L_0/dx;           % Lattice size, unit site

N=(1:mode)';
W=zeros(2,size(N,1));
W(1,freq)=a*sqrt(L/2);
k=(para(2)*para(1))./(2*(-4*para(1).*repmat(N,1,length(t)).^2*pi*pi/L/L));
lambda=-4*para(1).*N.^2*pi*pi/L/L;

a2=k.*(exp(2.*lambda*t)-1);
r2=(2*sum(a2,1)+(W(1,freq)*(exp(lambda(freq,1)*t)-1)/lambda(freq,1)).^2)/L;

end